Guidelines for Free-Energy Calculations Involving Charge Changes.

The Coulomb interactions in molecular simulations are inherently approximated due to the finite size of the molecular box sizes amenable to current-day compute power. Several methods exist for treating long-range electrostatic interactions, yet these approaches are subject to various finite-size-related artifacts. Lattice-sum methods are frequently used to approximate long-range interactions; however, these approaches also suffer from artifacts which become particularly pronounced for free-energy calculations that involve charge changes. The artifacts, however, also affect the sampling when plain simulations are performed, leading to a biased ensemble. Here, we investigate two previously described model systems to determine if artifacts continue to play a role when overall neutral boxes are considered, in the context of both free-energy calculations and sampling. We find that ensuring that no net-charge changes take place, while maintaining a neutral simulation box, may be sufficient provided that the simulation boxes are large enough. Addition of salt to the solution (when appropriate) can further alleviate the remaining artifacts in the sampling or the calculated free-energy differences. We provide practical guidelines to avoid finite-size artifacts.

Phase-Dependent Adsorption of Myelin Basic Protein to Phosphatidylcholine Lipid Bilayers.

The dense packing of opposite cytoplasmic surfaces of the lipid-enriched myelin membrane, responsible for the proper saltatory conduction of nerve impulses through axons, is ensured by the adhesive properties of myelin basic protein (MBP). Although preferentially interacting with negatively charged phosphatidylserine (PS) lipids, as an intrinsically disordered protein, it can easily adapt its shape to its immediate environment and thus adsorb to domains made of zwitterionic phosphatidylcholine (PC) lipids. As the molecular-level interaction pattern between MBP and PC lipid membranes suffers from scarce characterization, an experimental and computational study of multilamellar liposomes (MLVs) composed of 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) in the presence of bovine MBP is presented here. Calorimetric and temperature-dependent UV-Vis measurements identified DPPC pretransition temperature (Tp) and calorimetric enthalpy (ΔHcal) as the physicochemical parameters most responsive to the presence of MBP. Besides suggesting an increase in ß-sheet fractions of structured MBP segments as DPPC lipids undergo from the gel (20 °C) to the fluid (50 °C) phase, FTIR spectra unraveled the significant contribution of lysine (Lys) residues in the adsorption pattern, especially when DPPC is in the fluid (50 °C) phase. In addition to highlighting the importance of Lys residues in the MBP adsorption on DPPC lipid bilayer, employing salt bridges (SBs) and hydrogen bonds (HBs), MD data suggest the crucial importance of the orientation of MBP with respect to the surface of the DPPC lipid bilayer.

Deciphering the origin of the melting profile of unilamellar phosphatidylcholine liposomes by measuring the turbidity of its suspensions.

The elucidation of the thermal properties of phosphatidylcholine liposomes is often based on the analysis of the thermal capacity profiles of multilamellar liposomes (MLV), which may qualitatively disagree with those of unilamellar liposomes (LUV). Experiments and interpretation of LUV liposomes is further complicated by aggregation and lamellarization of lipid bilayers in a short time period, which makes it almost impossible to distinguish the signatures of the two types of bilayers. To characterize independently MLV and LUV of 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC), the latter were prepared with the addition of small amounts of 1,2-dipalmitoyl-sn-glycero-3-phosphatidylglycerol (DPPG) which, due to the sterical hindrance and negative charge at a given pH value, cause LUV repellence and contribute to their stability. Differential scanning calorimetry curves and temperature-dependent UV/Vis spectra of the prepared MLV and LUV were measured. Multivariate analysis of spectrophotometric data determined the phase transition temperatures (pretransition at Tp and the main phase transition at Tm), and based on the changes in turbidities, the thickness of the lipid bilayer in LUV was determined. The obtained data suggested that the curvature change is a key distinguishing factor in MLV and LUV heat capacity profiles. By combining the experimental results and those obtained by MD simulations, the interfacial water layer was characterized and its contribution to the thermal properties of LUV was discussed.

Perturbation Free-Energy Toolkit: An Automated Alchemical Topology Builder.

Free-energy calculations play an important role in the application of computational chemistry to a range of fields, including protein biochemistry, rational drug design, or materials science. Importantly, the free-energy difference is directly related to experimentally measurable quantities such as partition and adsorption coefficients, water activity, and binding affinities. Among several techniques aimed at predicting free-energy differences, perturbation approaches, involving the alchemical transformation of one molecule into another through intermediate states, stand out as rigorous methods based on statistical mechanics. However, despite the importance of free-energy calculations, the applicability of the perturbation approaches is still largely impeded by a number of challenges, including the definition of the perturbation path, i.e., alchemical changes leading to the transformation of one molecule to the other. To address this, an automatic perturbation topology builder based on a graph-matching algorithm is developed, which can identify the maximum common substructure (MCS) of two or multiple molecules and provide the perturbation topologies suitable for free-energy calculations using the GROMOS and the GROMACS simulation packages. Various MCS search options are presented leading to alternative definitions of the perturbation pathway. Moreover, perturbation topologies generated using the default multistate MCS search are used to calculate the changes in free energy between lysine and its two post-translational modifications, 3-methyllysine and acetyllysine. The pairwise free-energy calculations performed on this test system led to a cycle closure of 0.5 ± 0.3 and 0.2 ± 0.2 kJ mol-1, with GROMOS and GROMACS simulation packages, respectively. The same relative free energies between the three states are obtained by employing the enveloping distribution sampling (EDS) approach when compared to the pairwise perturbations. Importantly, this toolkit is made available online as an open-source Python package (https://github.com/drazen-petrov/SMArt).

Exploring the structure and dynamics of proteins in soil organic matter.

Alongside inorganic materials, water, and air, soil organic matter (SOM) is one of the major components of soil and has tremendous influence on the environment given its vital role in the carbon cycle. Many soil dwelling organisms like plants, fungi and bacteria excrete proteins, whose interaction with SOM is poorly understood on an atomistic level. In this study, molecular dynamics simulations were used to investigate selected proteins in soil models of different complexity from simple co-solvent molecules to Leonardite humic acids (LHA). We analyzed the proteins in terms of their structural stability, the nature and strength of the interactions with their surroundings, as well as their aggregation behavior. Upon insertion of proteins in complex SOM models, their structural stability decreased, although no unfolding or disruption of secondary structure was observed. The interactions of proteins and SOM were primarily governed by electrostatic forces, often in form of hydrogen bonds. However, also weaker van der Waals forces made a significant contribution to the total interaction energies. Moreover, we showed that even though the molecular structure and size of SOM molecules varied, the functional groups of SOM ordered around the protein in a similar pattern. Finally, the number of aggregates formed by proteins and SOM molecules was shown to be primarily proportional to the size of the latter. Strikingly, for varying protein net charges no changes in the formation of aggregates with the strongly negatively charged LHA were observed.

Vienna soil organic matter modeler 2 (VSOMM2).

Soil Organic Matter (SOM) plays an important role in several biogeochemical processes by directly affecting the microbial activity, soil aggregation, plant growth and carbon storage. Despite of its importance, our understanding of its composition and structure is still incomplete. Several experiments using elemental analysis, nuclear magnetic resonance (NMR) and mass spectrometry (MS) shed light on the structure of organic matter. In this context, the Vienna Soil-Organic-Matter Modeler (https://somm.boku.ac.at/) is a website that generates condensed phase computer models of Soil-Organic-Matter (SOM). Most of the data comes from standardized samples by the International Humic Substances Association (IHSS), which uses a specific methodology to extract organic compounds from soil, called humic substances. We have improved the modeler by increasing the pool of elemental units that compose our SOM molecules called building blocks, and also by implementing a genetic algorithm that increases the chemical and geometric diversity of the models. This allowed us to create models using the IHSS data as well as different types of soil. The webserver uses as an input principally the elemental and organic composition and offers input files needed to run molecular dynamic (MD) simulations of solvated and neutralized SOM within the framework of the GROMOS 54A7 forcefield and the GROMOS and GROMACS simulation packages.

Optimization of Alchemical Pathways Using Extended Thermodynamic Integration.

Thermodynamic integration (TI) is a commonly used method to determine free-energy differences. One of its disadvantages is that many intermediate λ-states need to be sampled in order to be able to integrate accurately over ⟨∂H/∂λ⟩. Here, we use the recently introduced extended TI to study alternative parameterizations of H(λ) and its influence on the smoothness of the ⟨∂H/∂λ⟩ curves as well as the efficiency of the simulations. We find that the extended TI approach can be used to select curves of low curvature. An optimal parameterization is suggested for the calculation of hydration free energies. For calculations of relative binding free energies, we show that optimized parameterizations of the Hamiltonian in the unbound state also effectively lower the curvature in the bound state of the ligand.

Molecular modelling of sorption processes of a range of diverse small organic molecules in Leonardite humic acid.

Soil organic matter (SOM) is abundant in the environment and plays an important role in several biogeochemical processes, including microbial activity, soil aggregation, plant growth and carbon storage. One of its key functions is the retention and release of various chemical compounds, primarily governed by the sorption process, which strongly affects the environmental fate of nutrients and pollutants. Sorption largely depends on the composition of SOM, as well as its structure, dynamics and the thermodynamic conditions. Although several approaches are available, experimental characterization of sorption mechanisms is not easy. Computational models for predicting sorption coefficients often require a wealth of experimental data for training and are only applicable to compounds and conditions related to the training dataset. Here, we use molecular dynamics (MD) simulations to study the sorption of a range of small organic compounds. As a model SOM system we use the standard Leonardite humic acid (LHA) sample, which physicochemical properties have recently been characterized computationally in detail. This model allowed us to estimate sorption propensities of the systems at two different hydration levels (water activities close to 0 and 1), showing a remarkable correlation with experimental data. Importantly, this molecular modelling approach based on perturbation free-energy calculations is rigorously derived from statistical thermodynamics and requires no experimental sorption data for training. It is therefore in principle applicable to any SOM model or thermodynamic condition. Moreover, the power of MD simulations to provide high-resolution insight into atomistic and molecular interactions was employed to explore how sorbate molecules associate with the LHA matrix and which contacts they form. The heteroatoms of both sorbate and sorbent play an important role and water molecules are identified as further key players in facilitating the sorption process. HIGHLIGHTS: Modelling of the sorption processes in soil organic matter at atomistic level.Rigorous, physics-based approach applicable to a range of SOM systems and conditions.Remarkable level of matching with experimental data with additional insight into the molecular mechanism.Interactions between the sorbate and local environment strongly affects the sorption process.

Toward Automated Free Energy Calculation with Accelerated Enveloping Distribution Sampling (A-EDS).

Free-energy perturbation (FEP) methods are commonly used in drug design to calculate relative binding free energies of different ligands to a common host protein. Alchemical ligand transformations are usually performed in multiple steps which need to be chosen carefully to ensure sufficient phase-space overlap between neighboring states. With one-step or single-step FEP techniques, a single reference state is designed that samples phase-space not only representative of a full transformation but also ideally resembles multiple ligand end states and hence allows for efficient multistate perturbations. Enveloping distribution sampling (EDS) is one example for such a method in which the reference state is created by a mathematical combination of the different ligand end states based on solid statistical mechanics. We have recently proposed a novel approach to EDS which enables efficient barrier crossing between the different end states, termed accelerated EDS (A-EDS). In this work, we further simplify the parametrization of the A-EDS reference state and demonstrate the automated calculation of multiple free-energy differences between different ligands from a single simulation in three different well-described drug design model systems.

Enhancing the promiscuity of a member of the Caspase protease family by rational design.

The N-terminal cleavage of fusion tags to restore the native N-terminus of recombinant proteins is a challenging task and up to today, protocols need to be optimized for different proteins individually. Within this work, we present a novel protease that was designed in-silico to yield enhanced promiscuity toward different N-terminal amino acids. Two mutations in the active-site amino acids of human Caspase-2 were determined to increase the recognition of branched amino-acids, which show only poor binding capabilities in the unmutated protease. These mutations were determined by sequential and structural comparisons of Caspase-2 and Caspase-3 and their effect was additionally predicted using free-energy calculations. The two mutants proposed in the in-silico studies were expressed and in-vitro experiments confirmed the simulation results. Both mutants showed not only enhanced activities toward branched amino acids, but also smaller, unbranched amino acids. We believe that the created mutants constitute an important step toward generalized procedures to restore original N-termini of recombinant fusion proteins.

Correcting electrostatic artifacts due to net-charge changes in the calculation of ligand binding free energies.

Alchemically derived free energies are artifacted when the perturbed moiety has a nonzero net charge. The source of the artifacts lies in the effective treatment of the electrostatic interactions within and between the perturbed atoms and remaining (partial) charges in the simulated system. To treat the electrostatic interactions effectively, lattice-summation (LS) methods or cutoff schemes in combination with a reaction-field contribution are usually employed. Both methods render the charging component of the calculated free energies sensitive to essential parameters of the system like the cutoff radius or the box side lengths. Here, we discuss the results of three previously published studies of ligand binding. These studies presented estimates of binding free energies that were artifacted due to the charged nature of the ligands. We show that the size of the artifacts can be efficiently calculated and raw simulation data can be corrected. We compare the corrected results with experimental estimates and nonartifacted estimates from path-sampling methods. Although the employed correction scheme involves computationally demanding continuum-electrostatics calculations, we show that the correction estimate can be deduced from a small sample of configurations rather than from the entire ensemble. This observation makes the calculations of correction terms feasible for complex biological systems. To show the general applicability of the proposed procedure, we also present results where the correction scheme was used to correct independent free energies obtained from simulations employing a cutoff scheme or LS electrostatics. In this work, we give practical guidelines on how to apply the appropriate corrections easily.

Polarization Effects in Simulations of Kaolinite-Water Interfaces.

Computational models of clay minerals and their interactions with the surrounding medium are highly valuable to study adsorption processes at an atomistic resolution, which may be relevant in different areas such as chromatography, environmental chemistry, and so forth. In this work, we analyzed the effect of the treatment of long-range interactions on the polarization of kaolinite-water interfacial systems in terms of structural, electric and dynamic properties, and hydrogen bonds. When using conventional three-dimensional (3D) Ewald summation, water molecules were more structured on the alumina interface of the kaolinite compared to simulations, in which the periodicity perpendicular to the plane was effectively removed. By applying an external electric field to the latter simulations, we were able to reproduce results obtained with the 3D Ewald summation. We proved that the induced polarization promotes adsorption processes of polar and charged species from aqueous solutions using a glyphosate molecule (widely used herbicide) as a particular example. The polarization phenomena raised from the simulation setup should be accounted for carefully in simulations of adsorption processes, which involve periodic images of polar interfaces. An inappropriate treatment of long-range electrostatic interactions can easily lead to artifacts and/or erroneous results.

The two cathepsin B-like proteases of Arabidopsis thaliana are closely related enzymes with discrete endopeptidase and carboxydipeptidase activities.

The genome of the model plant Arabidopsis thaliana encodes three paralogues of the papain-like cysteine proteinase cathepsin B (AtCathB1, AtCathB2 and AtCathB3), whose individual functions are still largely unknown. Here we show that a mutated splice site causes severe truncations of the AtCathB1 polypeptide, rendering it catalytically incompetent. By contrast, AtCathB2 and AtCathB3 are effective proteases which display comparable hydrolytic properties and share most of their substrate specificities. Site-directed mutagenesis experiments demonstrated that a single amino acid substitution (Gly336→Glu) is sufficient to confer AtCathB2 with the capacity to tolerate arginine in its specificity-determining S2 subsite, which is otherwise a hallmark of AtCathB3-mediated cleavages. A degradomics approach utilizing proteome-derived peptide libraries revealed that both enzymes are capable of acting as endopeptidases and exopeptidases, releasing dipeptides from the C-termini of substrates. Mutation of the carboxydipeptidase determinant His207 also affected the activity of AtCathB2 towards non-exopeptidase substrates, highlighting mechanistic differences between plant and human cathepsin B. This was also noted in molecular modeling studies which indicate that the occluding loop defining the dual enzymatic character of cathepsin B does not obstruct the active-site cleft of AtCathB2 to the same extent as in its mammalian orthologues.

Ion-induced modification of the sucrose network and its impact on melting of freeze-dried liposomes. DSC and molecular dynamics study.

Disaccharides play an important role in survival of anhydrobiotic organisms during extreme environmental conditions. A key protection feature is their capability to form the hydrogen bond (HB) network in a similar fashion as the one made by water. Since various ions also affect the HB network in completely hydrated systems, it is of a great interest to understand how they impact preservation when incorporated in a disaccharide network. To address this, we employ a combination of experimental and modeling techniques to study behavior of multilamellar 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) liposomes freeze-dried with sucrose in presence of NaCl or NaH2PO4·H2O at various concentrations (0.01-1M). Differential scanning calorimetry (DSC) was employed in order to determine the cooperative unit size (CUS), the number of lipid molecules that constitute a domain of cooperative motion in the liposome, and the melting temperature (Tm). In the absence of salt CUS was estimated to be 122±12, whereas in the presence of NaCl CUS increases more (347±34 for c=1M) than for NaH2PO4·H2O (193±26 for 1M). When it comes to Tm, the situation is reversed; NaCl induces increase by about 1K, while NaH2PO4·H2O by about 10K. These findings clearly demonstrate how different interaction forces-hydrogen bonding, charge pairing, and van der Waals interactions between acyl chains-affect CUS and Tm. Their interplay and contribution of particular interaction was further analyzed with molecular dynamics (MD) simulations. This analysis demonstrated that the HB network of DMPC and sucrose is partially disrupted in the presence of NaCl ions, and even to a greater extent in the case of NaH2PO4·H2O ions. Notably, H2PO4- ions outcompete and replace the sucrose molecules at the DMPC surface, which in turn alters the nature of the DMPC-surrounding interactions, from a weaker HB-dominated to a stronger CP-dominated interaction network.

Molecular Dynamics Simulations of the Standard Leonardite Humic Acid: Microscopic Analysis of the Structure and Dynamics.

Humic substances (HS) are abundant in the environment and play an important role in a number of biogeochemical processes including microbial activity, soil aggregation, plant growth, the retention and release of nutrients, the environmental fate of pollutants, and carbon storage. They are flexible, relatively small molecules forming supramolecular structures through weak interactions. Despite the great importance of understanding their behavior at the atomic level, computational modeling, a premier high-resolution technique providing great level of detail, has been surprisingly little-employed to study humic substances. Here, we use the recently developed Vienna Soil Organic-Matter Modeler to create representative models of a real HS sample, the standard Leonardite humic acid. Molecular dynamics simulations were used to probe the structure and dynamics of the system at a range of hydration levels. The studied systems were characterized in terms of their physicochemical properties, including density, dielectric properties, hydrogen bonding, etc. Moreover, the strength of sorption was estimated for three small organic compounds: benzaldehyde, propan-2-ol, and acetone. Strikingly, the HS models were validated against experimental data showing a remarkable agreement with calculated properties. Finally, we make the equilibrated models of the standard Leonardite humic acid, together with corresponding force-field parameters, available at the Vienna Soil Organic-Matter Modeler.

Effect of Oxidative Damage on the Stability and Dimerization of Superoxide Dismutase 1.

During their life cycle, proteins are subject to different modifications involving reactive oxygen species. Such oxidative damage to proteins may lead to the formation of insoluble aggregates and cytotoxicity and is associated with age-related disorders including neurodegenerative diseases, cancer, and diabetes. Superoxide dismutase 1 (SOD1), a key antioxidant enzyme in human cells, is particularly susceptible to such modifications. Moreover, this homodimeric metalloenzyme has been directly linked to both familial and sporadic amyotrophic lateral sclerosis (ALS), a devastating, late-onset motor neuronal disease, with more than 150 ALS-related mutations in the SOD1 gene. Importantly, oxidatively damaged SOD1 aggregates have been observed in both familial and sporadic forms of the disease. However, the molecular mechanisms as well as potential implications of oxidative stress in SOD1-induced cytotoxicity remain elusive. In this study, we examine the effects of oxidative modification on SOD1 monomer and homodimer stability, the key molecular properties related to SOD1 aggregation. We use molecular dynamics simulations in combination with thermodynamic integration to study microscopic-level site-specific effects of oxidative "mutations" at the dimer interface, including lysine, arginine, proline and threonine carbonylation, and cysteine oxidation. Our results show that oxidative damage of even single residues at the interface may drastically destabilize the SOD1 homodimer, with several modifications exhibiting a comparable effect to that of the most drastic ALS-causing mutations known. Additionally, we show that the SOD1 monomer stability decreases upon oxidative stress, which may lead to partial local unfolding and consequently to increased aggregation propensity. Importantly, these results suggest that oxidative stress may play a key role in development of ALS, with the mutations in the SOD1 gene being an additional factor.

Are current atomistic force fields accurate enough to study proteins in crowded environments?

The high concentration of macromolecules in the crowded cellular interior influences different thermodynamic and kinetic properties of proteins, including their structural stabilities, intermolecular binding affinities and enzymatic rates. Moreover, various structural biology methods, such as NMR or different spectroscopies, typically involve samples with relatively high protein concentration. Due to large sampling requirements, however, the accuracy of classical molecular dynamics (MD) simulations in capturing protein behavior at high concentration still remains largely untested. Here, we use explicit-solvent MD simulations and a total of 6.4 µs of simulated time to study wild-type (folded) and oxidatively damaged (unfolded) forms of villin headpiece at 6 mM and 9.2 mM protein concentration. We first perform an exhaustive set of simulations with multiple protein molecules in the simulation box using GROMOS 45a3 and 54a7 force fields together with different types of electrostatics treatment and solution ionic strengths. Surprisingly, the two villin headpiece variants exhibit similar aggregation behavior, despite the fact that their estimated aggregation propensities markedly differ. Importantly, regardless of the simulation protocol applied, wild-type villin headpiece consistently aggregates even under conditions at which it is experimentally known to be soluble. We demonstrate that aggregation is accompanied by a large decrease in the total potential energy, with not only hydrophobic, but also polar residues and backbone contributing substantially. The same effect is directly observed for two other major atomistic force fields (AMBER99SB-ILDN and CHARMM22-CMAP) as well as indirectly shown for additional two (AMBER94, OPLS-AAL), and is possibly due to a general overestimation of the potential energy of protein-protein interactions at the expense of water-water and water-protein interactions. Overall, our results suggest that current MD force fields may distort the picture of protein behavior in biologically relevant crowded environments.

Trichoplaxin - a new membrane-active antimicrobial peptide from placozoan cDNA.

A method based on the use of signal peptide sequences from antimicrobial peptide (AMP) precursors was used to mine a placozoa expressed sequence tag database and identified a potential antimicrobial peptide from Trichoplax adhaerens. This peptide, with predicted sequence FFGRLKSVWSAVKHGWKAAKSR is the first AMP from a placozoan species, and was named trichoplaxin. It was chemically synthesized and its structural properties, biological activities and membrane selectivity were investigated. It adopts an α-helical structure in contact with membrane-like environments and is active against both Gram-negative and Gram-positive bacterial species (including MRSA), as well as yeasts from the Candida genus. The cytotoxic activity, as assessed by the haemolytic activity against rat erythrocytes, U937 cell permeabilization to propidium iodide and MCF7 cell mitochondrial activity, is significantly lower than the antimicrobial activity. In tests with membrane models, trichoplaxin shows high affinity for anionic prokaryote-like membranes with good fit in kinetic studies. Conversely, there is a low affinity for neutral eukaryote-like membranes and absence of a dose dependent response. With high selectivity for bacterial cells and no homologous sequence in the UniProt, trichoplaxin is a new potential lead compound for development of broad-spectrum antibacterial drugs.

A systematic framework for molecular dynamics simulations of protein post-translational modifications.

By directly affecting structure, dynamics and interaction networks of their targets, post-translational modifications (PTMs) of proteins play a key role in different cellular processes ranging from enzymatic activation to regulation of signal transduction to cell-cycle control. Despite the great importance of understanding how PTMs affect proteins at the atomistic level, a systematic framework for treating post-translationally modified amino acids by molecular dynamics (MD) simulations, a premier high-resolution computational biology tool, has never been developed. Here, we report and validate force field parameters (GROMOS 45a3 and 54a7) required to run and analyze MD simulations of more than 250 different types of enzymatic and non-enzymatic PTMs. The newly developed GROMOS 54a7 parameters in particular exhibit near chemical accuracy in matching experimentally measured hydration free energies (RMSE=4.2 kJ/mol over the validation set). Using this tool, we quantitatively show that the majority of PTMs greatly alter the hydrophobicity and other physico-chemical properties of target amino acids, with the extent of change in many cases being comparable to the complete range spanned by native amino acids.

Vienna-PTM web server: a toolkit for MD simulations of protein post-translational modifications.

Post-translational modifications (PTMs) play a key role in numerous cellular processes by directly affecting structure, dynamics and interaction networks of target proteins. Despite their importance, our understanding of protein PTMs at the atomistic level is still largely incomplete. Molecular dynamics (MD) simulations, which provide high-resolution insight into biomolecular function and underlying mechanisms, are in principle ideally suited to tackle this problem. However, because of the challenges associated with the development of novel MD parameters and a general lack of suitable computational tools for incorporating PTMs in target protein structures, MD simulations of post-translationally modified proteins have historically lagged significantly behind the studies of unmodified proteins. Here, we present Vienna-PTM web server (http://vienna-ptm.univie.ac.at), a platform for automated introduction of PTMs of choice to protein 3D structures (PDB files) in a user-friendly visual environment. With 256 different enzymatic and non-enzymatic PTMs available, the server performs geometrically realistic introduction of modifications at sites of interests, as well as subsequent energy minimization. Finally, the server makes available force field parameters and input files needed to run MD simulations of modified proteins within the framework of the widely used GROMOS 54A7 and 45A3 force fields and GROMACS simulation package.